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1.0 SDMfART AND INTRODUCTION 


Although great strides have been made in the area of computational fluid 
mechanics in the past 40 years, it is probably fair to say that the progress 
has principally been due to enormous strides in the capacities of machinery 
rather than in the science of numerical methods. Problems of stability, 

conservation, and the treatment of corner flows inviscidly combined to 
restrict the ability of numericists to solve many problems of current 
interest. 

In a previous discussion pertaining to one-dimensional compressible 

f 1 ^ 

unsteady flow by Prozan'- , a variational principle for fluid mechanics was 
hypothesized. The excellent stability characteristics and absolutely 
conservative properties of that solution encouraged the extension of the 
investigation to multi-dimensional flow. The discussion in this report only 
deals with two-dimensional flow for square elements, but it is felt that the 
extension to three-dimensional flow and/or general element shapes presents no 
fundamental problems. 

Introduction of the variational principle utilized herein indicates that 
the satisfaction of equations of motion which have been heretofore considered 
necessary and sufficient have only been necessary but not sufficient to 
guarantee a successful solution. Although successful solutions have been 
produced, they have been produced at the expense of conservation. Because (at 
least for explicit methods) absolutely conservative methods exhibit 
undesirable behavior with regards to stability, the conservative properties 
have been traded off against stability. 

The key features of the methodology presented in the ensuing discussion 

are: 

o absolutely conservative 

o no explicit additions of artificial viscosity 

o absolutely stable - stable time step internally 
controlled - operates near the CFL 


o 


dynamically differences based on local conditions 

o can be run with sparse grids for engineering or 
parametric solutions - can also be run with fine grid 
for more precise answers 

o runs first time every time - no operator interaction 
nor special expertise is required. 

In short, the methodology may be used to provide a quick response engineering 
analysis to whatever level of accuracy is required. 
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2.0 DISCUSSION 



2.1 THE VARIATIONAL PRINCIPLES 

For an unreacting, inviscid fluid, the Second Law of Thermodynamics 
states that 


+ 7 • pqs)dV > 0 (1) 

where p is the density, s is the intensive entropy, q is the velocity vector, 
and V is the volume of the domain. In the following discussion we will assume 
that the equations of motion are satisfied. Now equation (1) states that the 
internal entropy (extensive) production is greater than the net flux for all 
time. Since the production exceeds the flux then the accumulation of entropy 
must constantly increase. If we define the outer boundaries of the problem 
far enough away from the area of interest (or consider a closed system) then 
we can consider the inlet/outlet flux terms as constant. Rewriting 
equation ( 1 ) ; 


Jv ^ dv - pqs . dA + pqs . dA > 0 (2) 
enables us to see that for an arbitrarily large time interval we have; 

dV + C > 0 (3) 
where C is a constant (and is zero for a closed system). 

It was previously stated that the entropy contained within the system 
must continually increase and since it cannot increase without limit then it 
must approach a maximum. Furthermore, since all unsteady fluctuations produce 
an entropy increase, however small, then the maximum entropy must occur at 
steady state. 

If we view the numerical approximation as a process acting on the system 
we may conclude that the inexactitude of the approximation must be such that 
the predicted entropy is greater than or equal to the exact solution entropy 
in order to guarantee stability. 
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To clarify this conclusion we must state the equations of motion; 


3t 3x 



0 


( 4 ) 


where 

e = (p, pu, pv, pE); 
f = (pu, p+pu^, puv, puH); 
g = (Pv, Puv, p+Pv^, pvH); 

and D represents the continuity, axial momentum, lateral momentum, and energy 
equations, respectively, and where u is the axial velocity and v the lateral 
velocity. 


Equation ( 1 ) may now be expanded as 


Q^D. dV > 0 (i=1,2,3,4) 


(5) 


where 


Q. = 1^ 

1 3e. 

1 

Obviously, if the equations of motion are satisified everywhere within the 
system (D^ = 0), then equation (5) is satisfied. Under the numerical 
approximation, however, the equations of motion are not satisfied exactly but 
only in some approximate sense. Equation (5) may or may not be satisfied, 
therefore, even if the equations of motion are satisfied in the approximate 
sense . 

Therefore, formulations which may be entirely acceptable in terms of the 
conservation equations may be unacceptable to (not satisfy) equation (5). It 
becomes clear then that satisfaction of the equations of motion is a necessary 
but not sufficient condition to satisfying the physical principles governing 
the motion. 

It was earlier stated that the equations of motion would be assumed 
satisfied for discussion purposes. The approach taken to accomplish this will 
now be discussed. 
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The equations of motion act as equality constraints on the formation of 
extensive system entropy. Using the method of Lagrange to state the 
constrained functional yields; 

S = (^ + V • pis3dV + X.D.dV > 0 (6) 

where S is the constrained extensive system entropy production and are the 
Lagrange multipliers. 

Differentiation 
to 9e^/9t yields; 


of the constrained functional with respect 


-i- = fv 


( 7 ) 


so that if the equality is chosen (exact solution) in equation (6) then we may 
conclude that the 


X. = - Q. (8) 
1 1 

This, of course, is consistent with equation (5) which states that the 
constraints degenerate in the limit and that the Euler-Lagrange equations for 
the functional described in equation (1) are the equations of motion. 

There could be, however, a whole class of functions which satisfy the 
above criterion (although only the one described above is embodied in a known 
principle of physics). In addition, we have no guarantee that a tractable 
solution can be devised, based on the principle, which offers a substantial 
advantage over current better known approaches. Of course, the success of the 
reference ( 1 ) solution for one-dimensional flow was an encouraging first step 
in the investigation. 

Before proceeding with the numerical procedure which was developed to 
ascertain the applicability of the stated variational principle to numerical 
solutions in multi-dimensional flow, an interpretation of the Lagrange 
multipliers appearing in equation (6) is in order. 
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The integral constraint term of equation (6) looks exactly like a Method 
of Weighted Residuals statement and in point of fact enters into the problem 
formulation in an identical fashion as the quasi-variational definition of 
free parameters used by Oden^^^ for the purpose of assembling finite elements 
into the entire computational domain. 

Before concluding this section of the development it should be pointed 
out that for a closed system at constant temperature and pressure the 
principle stated does degenerate to Gibbs^^^ criteria for equilibrium of a 
system of known energy. This observation suggests that the variational 
statement applies also to reacting flows. It further lends credence to a 
hypothesis by Prozan^^^ which contends that there is a coupling between 
chemistry and flow that is not accounted for in the current usage of Gibbs 
free energy. 

2.2 THE METHODOLOGY 

2.2.1 Construction of the numerical model 

Proceeding now with the construction of the model, we must first define 
p = ^ {pE - -^(pu^ + pv^)}, H = E + ^ 

TT ^ 


and 


s = c^ (In p - Y P) where p = pRT 


Also , 


9ps 

9p 


(Ts - H + q^)/T; 


9ps 

3pu 


-u/T 


9ps 

9pv 


-v/T; 


8ps 
3pE = 


1/T 


If the domain of interest is subdivided into a number of elements then we may 
write 
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S = I L 


e=1 


V ' 3t 3x ^ 3y ^ 


E 3e. 3f. 3g. 

+ I/v + -§J + -§7) dV„ > 0 

e=1 e 


e - 


(9) 


y 


Sketch 1 - Element Description 


K 

— ► X 


where the corners of the element are nodes in the assembled system. Due to 
the exploratory nature of this development, simplifying assumptions were made 
which greatly reduced the development lead time but must be removed before a 
general numerical model can be devised. Thus the element was assumed square 


with sides of length one. Therefore we transform; 

x=Xj+^;y=yj+n (10) 

The following shape functions were also assumed; 

e = e^(l-C)d-n) + e^jCd-D) + + ejy(1-^)n (11a) 

f = fjd-Od-n) + fjj5d-n) + + fjy(i-^)n (iib) 

g = gjd-Od-n) + gjj5d-h) + gjjjSn + gjyd-^)n die) 

X = XjA^(5,n)+XjjAjj(5,Ti)+XjjjAjjj(^,Ti)+Xj^jy(^,Ti) did) 


9 
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Notice that the "weight" functions Aj, Ajj, Ajjj, Ajy have not been specified 
at this juncture. 

Substitution of the above relations into a single element of the 
constraint terms of equation (9) yields, in general, an implicit analog. It 
is of interest, however, to create an explicit analog. If we choose the 
weight functions orthogonal to the shape functions, an explicit analog 
results. 


h - - Tl) 

hi - 

^iii ■ '’ill (3 ■ 5) (3 - n] 

^xv ~ '’iv * ^^^"3 ** 

(This orthogonality observation is the basis for the reference (5) analysis 
where a more complete discussion may be found). This step ("explicitizing" 
the solution) introduces undetermined coefficients bj, bjj, bjjj, bjy which 
may be chosen in various ways and are therefore critical in controlling the 
solution behavior. 

Substitution of the shape and weight functions into a single element 
integral of the constraint term in equation (9) yields, after integration. 


(12a) 

(12b) 

(12c) 

(12d) 


%) ' ®I‘^^II“^I'^®IV"®I ^ 

^II^I:^^®II'^^II"^I'^®III"®II^ 
^III^II:^^®III'^^III"^IV'*‘®III"®II^ 
* ^IV^I1^^®IV'^^IIl”^IV‘^®IV“®I^ 


(13) 


The entire constraint term will now be formed. Consider a small subset 
of the assembled region. 


8 


7 8 9 


9 ^ 

© 

9 ...■■^ 

@ 

4 

© 

1 i 

5 6 

® 

t i 


1 2 3 


Sketch 2 - Assembly Description 


The constraint term integral becomes 


'b* b'" ••• " 

* ^2 - - 

* ^54 i ^^®5 * ^5 ■ ‘'4 ®5 

H. * fj - f2 ■. gj - 

* *’6 - ^5 * 86 

* ^5*’]0®5 * ^6 ■ ^5 * ®8 “ 

^ * ‘'9 - *"8 * 89 

* ^ 4'^]^^®4 + ^5 ■ *’4 87 - 

* ^8*’li:^^®8 * *'8 ■ *'7 83 


g^) + 

X,,b__ C©/, + 

2 2 

^2 - 

^1 " 

«5 - 

82) 

- 83) 

* ^ 4 *’l \^^®4 

- "5 

- ^ 

+ 64 

- 

g2> + 


^3 ■ 

^2 * 

S6 - 

g 3 ) 

- 83) 


* ^6 

- ^5 


- 83) 

85) 


^6 - 


«9 - 

85’ 

- 86> 

^ ®8 

"*■9 

- "8 

"■ «8 

- 85) 

84) + 


^5 - 

f 4 + 

®8 - 

8^ ) 

■ 


" *'8 

- ^7 

+ Sj 

- 84) 


s 0 


(14) 
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Notice that all terms with coefficients of are explicitly present in 
equation (14). The same may or may not be true for the other values of 
X-] - Xg depending on whether or not they are on boundaries. 

A digression from the element derivation is taken here to explicitly 
recreate the nodal analog. This step is not actually in the computation, but 
is undertaken at this point with the intention of clarifying the analog. 

Collecting coefficients of X^ yields, for the constraint term. 






■ ■ ^5 ®8 ■ 

- g^)} + ... = 0 



(15) 


In the quasi-variational formulation of Oden, X^ is considered an 
arbitrary parameter, therefore its coefficient must be zero. In this 
formulation, we recapture the constraint equation by differentiation with 
respect to the Lagrange multipliers. The result is the same for either 
interpretation; i.e. 


(b.___ + + b + b )e + b (f - f,, + g^. - g„) 




( 16 ) 


The above relation is the assembled nodal analog used for a general interior 
point. In the actual procedure this equation is not formulated and coded but 
rather is constructed as an assembly. Equation (16) is presented here, 
however, in its assembled form so that its features may be examined. 


Now at this point all of the b factors appearing in equation (16) are 
arbitrary. Choosing all parameters equal to 1/4 creates a centered space 
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backward time analog. Other choices give other analogs (perhaps not all of 
which have been attempted by previous researchers). 

Now that the analog has been explicitly illustrated, a return will be 
made to the element derivation. 


Collecting all coefficients of f^ 
volume integral results in 


JV 3t IV(g 


+ (b 


+ b^ + b. 




and g^ which will appear in a system 


+ . . . — 0 


( 17 ) 


If the conservation equations are to apply over the entire domain, then 
the contribution of f^ and g^ must disappear. Therefore, for integral 
conservation , we have ; 


° 

° 


(I8a) 

(I8b) 


To preserve the volume of the system we must have (contributions of 
subvolume pieces add up to the system volume); 


bj + b^j 




IV 


= 1 


(18c) 


It is desirable to have a scheme which defines the b parameters on an 
element level without regard to the surrounding elements. At least one such 
definition is, for each element, 


II 


1 

2 * 


*’lll * *’iv 


I . 

2 ’ 


‘'l = ‘’ill ' 


•’ll = *’IV 


( 19 ) 


which leaves us with one free parameter per element. This selection is 
modified for boundary elements in such a fashion that all conservation 
requirements are met while exactly satisfying boundary conditions. 


The definition of parameters found in equation (19) solves a wide range 
of very difficult problems but did exhibit instabilities at very high shock 
tube and diaphragm burst pressure ratios (1000/1). In view of the fact that 
the one-dimensional solution did successfully handle even higher pressure 
ratios (greater than 1000/1) it seemed that a superior technique could be 
created. Fortunately this was the case. 

If the original differential equation is split 



and the identical derivation is followed, we find that the coefficients of f^ 
and g^ are; 

'■to* '“to* 'to* '“to* 'to* ''to* °“to’ " 

where equation (21a) is identical to equation (l8a) but with a substitution of 
a for b while equation (21b) is identical to equation (I8b) with a 
substitution of 8 for b. 

The system volume integral is preserved if 

OLj + otjj + + Ojy = 1 (22a) 

3j + + 3 jy = T (22b) 

Equations (21) and (22) are satisfied if 

“l “ll ■ 2 ' °^II " 2 ’ ^I ^IV ” 2 ' ^II ^III ■ 2 

Again, only for discussion purposes, we write the nodal analog. 
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(eAV)^ 


where 




( 24 ) 


<•“>5 ■ <”lHj,* «Ife- «^- ”l^«*x>5 

* "“to* *”fi)* ®to* "‘to’"'’' "” 

Comparison of equations (24) and (25) with equation ( 16 ) shows that we have 
sacrificed nothing in the substitution of a and 3 for b. As we shall see 

later a great deal has been gained since the scheme based on a and 3 is much 
more powerful. 


Inspection of equations (22) and (23) reveals that there are now four 
free parameters (instead of only one using equations (19)). If we define 


II 

""s ’ 

II 

H 

¥ (1 - 

(26a) 

“iv = ¥ 

% ’ 

^‘lll 

II 

1 

(26b) 

h-i 

’ 

II 

> 

H 

OQ 

5 - V 

(26c) 

’ 2 

3g ; 


= 2 (1 - Bg) 

(26d) 


where otg, Ojg, 3yj, 3g may take on values between 0 and 1 (values less than zero 
or greater than one will generate negative volumes and are therefore avoided). 


It is necessary at this time to introduce boundary conditions. The 
computer code which was produced under this effort utilized a scheme which 
will be explained here. Although the identification scheme is not central to 
the theory it will simplify the explanation of the necessary treatment. 
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At each node (corner of a square element) in the system, a boundary 
condition code word is input. This code word consists of four digits (0 or 1) 
where each digit corresponds to an independent variable. The variables are 
taken in the order P, Pu, Pv, and PE. Zero means the parameter cannot change 
while one means that it is free. Thus 1011 identifies a tangency condition 
where pu cannot change but pv can while 1101 is an axial slip condition. A 
fixed inlet then is 0000 . An interior node is 1111. Consider then the 
element shown below. 


3 Sketch 3 - Free Parameters 

W 


a 

N 



Suppose now that the property at nodes I and II cannot change, i.e., 
eAv = 0. Unless this is a corner element then the element to the right and 
to the left will have the same condition. Further, let us assume that corner 
I is node 5 in equations (24,25). Obviously, regions © and (J) do not exist so 
that under those conditions equations (24,25) become. 


(eAV)^ 


where 




+ a 


T 


(6 


b 


)(gg - gg)} (27a) 


(eAV)^ = (a^ + 


b 


^ b 


®TT 


Apparently 


(27b) 


a = a =3, = =0 

for if they did not, (eAV)^ would have a value and (e) would change. We can 

D 

see that, within the element. 
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a, = = 0 ; 


«III «iv = ^ 



II 


= 0 


HI 


= a 



while maintaining absolute conservation , that there is only one free parameter 

(20j^). 


There are many other situations which can occur for corner nodes, nodes 
on left and right-hand walls. These situations may be derived by the 
interested reader but a treatment of each situation would considerably 
lengthen this discussion and is therefore omitted. The final scheme which 
treats all conditions is stated later. 

It is now possible to discuss the selection of the values for ag, 
and 6 g, which control the stability of the solution. Since the constraint 
equations have been satisfied we may write, for the system entropy production, 


S = = n«i.(eAV)j . 

e =1 e e =1 1 i i i 


^ QlII.(eAV)m + } 

1 1 1 L e 

where there is a contribution from each equation (i). Examining the 
contribution from a single equation within a typical element; 




or 




BgCBii - Bjjj) + 2 (Ajj + Ajjj + + Bj^j) 
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where 



= - : 


- Qj(gjv ~ Sj) 

(28a) 

Aji = 

^II^^II ” ^I^ ’ 

®ii = 

" *^11^^111 " ®ll^ 

(28b) 

Ajj-i = - 


®iii 

= - - gjj) 

(28c) 

A 

IV " 

" ^IV^^III ~ ^IV^ 

’ ®IV 


(28d) 


as. 

; 

aOg I II ’ 


^IV “ ^III 

(29a, b) 


as. 

3^6^ = - ®IV = 

93g 

®II " ®III 

(29c, d) 

The above equations 

(29) tell us how to 

change 

the free parameters 

in order to 

promote entropy production but give 

us no 

information on what 

value the 


parameters should have. Since the range of parameters is (0,1) we choose the 
parameters such that the most entropy is produced. Therefore, assuming that 
ctg, Oj^j, 3^^, Bg set to zero to begin with, then 


- 

° “s = * 

(30a) 

^IV “ 

'‘ill > 0 «N = 1 

(30b) 

"l - 

®IV ^ ° ' 

(30c) 

®II 

Bjii > 0 ^ gg = 1 

(30d) 


Notice that there is no guarantee that the system entropy production will 
exceed the flux but that within the constraints of the equations of motion, 
boundary conditions, and element interactions that the production rate is 
maximized . 

The corner node a, 3 then computed from equations (26). The boundary 
condition digits for this particular equation are summed. If the sum equals 
four then all nodes are free. If the sum is two or one, an alteration of the 
a,B is required to satisfy boundary conditions. 
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If the sum is two, then 


“j - 


6ji - ^ 0 


(31a) 


a. ^ 6. a. 
J 3 3 


*II - «I = 0 


(31b) 


while if the sura is one; 


“j " 


(31c) 


where 6 is the (0,1) digit signifying free or fixed and j is the node 
(I, II, III, IV). 

To asserable the element into the system the following procedure is used. 


(cAV)^ = (eAV)^ - 




(eAV)^ = (eAV)^ - ^II^^IIl“^IV^ 

3iii(giii-gll) 


4 IV III 3 


1 I 


II 2 


(eAV)' = (eAV)^ - Oj.{fjj-fj) 


^I^®IV"%^ 


(eAV>2 = (eAV)^ - 

- 3j^(gjjj-gjj) 


Sketch 4 - Derivative Assembly 
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In the diagram above, the nodes in the full domain numbering system are 1, 2, 
3, and At the beginning of each time step, the (eAV) for each node are set 
to zero and a sweep through each element is made adding in contributions from 
each corner of the element. In this fashion, the explicit analog is 
constructed dynamically at each time step. 

2.2.2 Expansion Corner Treatment 

Another very important facet of the development was the expansion corner 
treatment. Due to the very simple nature of the element shape and time and 
core restraints on the machine, the following corner procedure was developed. 

The element analog in the expansion corner region was assembled (on 
paper) and examined. Having only one node at a corner which was subject to 
axial slip on one side and lateral slip on the other, the only acceptable 
conclusion to satisfy both conditions was that the corner is a stagnation 
point. While this is certainly true in the real world, on the scale of the 
element sizes and nodal densities which we must accept in the development, a 
stagnation treatment was out of the question. A scheme was devised whereby 
special boundary conditions were enforced which allowed the corner to be 
treated as axial slip on one side, free at the corner, and lateral slip on the 
other. To clarify this, consider the assembly belov;. 


8 9 



Sketch 5 - Expansion Corner 
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Node 4 has boundary condition code of 1101 while node 8 is 1011 and node 5 is 
1111. A special condition is introduced where node III of element (J) is 1101 
and node I of element is 1011. Now if 6^ and 6^ are normally taken as unity 

for any node but are reset to the second and third digits of the special 
condition at node III of element (l) and node I of element Q) then the corner 
satisfies slip on both sides and is absolutely conservative as long as we 
redefine 


2 

f = 6pu, p + 6pu , 6 puv, 6 puH 
u u ’ uv ’ u 

2 

g = 6^pv, p + 6^pv , 6^pvH 

where = 0 at the corner (common) node. 

Since the pressure is continuous at the corner, the treatment leaves 
something to be desired for supersonic flow but seems to be an excellent 
solution for subsonic flow. 

To this point the solution has been first order time. A second order (or 
higher) method may be readily evolved. The statements in equations (1), (2) 
and (6) may also be made over an arbitrary time interval. For instance, 
equation (6) becomes 


(32a) 

(32b) 


S(t2) - S(ti) = S dt 
while our differential equation becomes 



The second order time solution was effected by letting; 




8 = - C2g^_^ 


(33) 


(34) 


(35a) 

(35b) 


where 
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and understanding 

interval At, , 
k 

To advance the assembled system nodal quantities one time step we let 

®k+1 (®^V)Atj^/Av (36) 

where Av is the assembled system volume associated with that node. To arrive 
at this system volume we associated 1/4 of each element surrounding the node 
with that node. A concave corner node for example has Av = 1/4 while a wall 
node has Av = 1/2 and an expansion corner has Av = 3/4 and finally an interior 
node has Av = 1 . 

2.3 the computer program 

A computer program was written to test the new procedure on an 

Alpha -Micro 100 micro computer at Continuum. Due to severe storage 
limitations, the largest grid attainable was 13 x 13 nodes (12 x 12 
elements). The computational sequence is 

1. Given the current state of the system (P, Pu, Pv, PE are 

known, at all nodes) compute f, g, Q at each node. 

Null eAv for all nodes. 

2. Sweep through all elements performing the following steps 
for each element and each equation for that element 

o Compute Aj,...,Ajy, using (28) 

o Compute oi^, 3^, 3g using (30) 

o Compute using (26) 

o Enforce boundary conditions using (31) 

o Assemble contribution into nodal system as in 
Sketch 4. 

3* Step forward in time using (36). 


, At, 
« . 1 k 

^1 = ^ 2 At 


k-1 


P I k 

^2 " 2 At 


k-1 


that (eAv) is the average change over the time 
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4 . Return to Step 1 . 

The CFL condition was used to compute the time step. An input quantity was 
used to define the maximum fraction of the CFL at which the solution was 
allowed to operate. The inequality constraint of non-zero and non-negative 
pressure and density was also enforced. If such a situation occurred, the 
time step was halved. For each successful time step, the time step was 
increased by M per cent (doubles every 20 steps) up to the limiting step. 

Each node was supplied with a boundary condition code. This information 
was principally used to give tangency conditions but also was used on density 
and energy if it was desired , for example , to keep an inlet fixed . 

An option was provided for first or second order time expansion. All 
results presented used second order time, but first order time also apparently 
works very well (this statement is qualified since not very many cases were 
run first order). 
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3.0 C»MPOTATIONAL RESULTS 


The first tests of the new computer program were one-dimensional shock 
tube and diaphragm burst problems. Due to the splitting technique outlined 
previously, the method degenerates exactly to one-dimensional flow formulation 
for these cases. Absolute conservation was achieved. The same results 
occurred whether the problems were run left to right, right to left, top to 
bottom, or bottom to top. No results are shown here since they were 
previously documented in reference ( 1) . 

Figure 1 illustrates the results of a problem which was run to determine 
corner conservation and long-term system entropy production. Absolute 
conservation was noted. The system extensive entropy did not monotonically 
increase but rather increased rapidly for the first 50 time steps and then 
oscillated with a general upward trend for the remainder of the problem. 

Figure 2 shows the velocity vectors resulting from the impulsive start of 
a Mach 1.2 jet (y = 1.4) exhausting into a quiescent background. The solution 
was run for 800 time steps at which time conditions were changing very 
slowly. The outflow boundaries of the problem were free. No instability 
occurred. 

Figure 3 illustrates flow over a rear-facing step and the results are 
shown at 200 time steps. Once again we have a Mach 1.2 jet (y = 1.4). This 
run was not an evolutionary run since the entire field was set to uniform 
parallel flow except at the back side of the step where the lateral slip 
condition was enforced. An inviscid recirculation cell has occurred along the 
back side of the step. No instabilities were noted (at least to this stage of 
the solution). 

Figure 4 gives the results for a Mach 1.2 jet (y r 1.4) into which a 
small high momentum jet is being injected. The Mach number of the secondary 
jet is 2.4 (y = 1.4). The solution flow vectors may be found in the figure as 
well as pressure distributions on the upper and lower walls. A high pressure 
region occurs just upstream of the secondary jet and a very low pressure 
region occurs downstream. The high pressure region along the upper wall has a 
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peak shifted downstream of the secondary jet position. The solution seemed 
very well behaved and was changing only slowly at 100 time steps. 

In general it can be stated that the method has demonstrated the 
following features: 

o No instabilities have occurred, 

o Absolute conservation, no matter what the complexity. 

o Successfully handles implusive starts of 1000/1 
pressure ratios for diaphragm burst and imposed jets. 

o Stability for closed systems approaching or achieving 
steady state. 

Due to the extremely long running time on this small machine, we cannot 
absolutely conclude that the solution is stable at very long run times when 
the solution approaches steady state. Problems of this nature must be 

investigated on a more powerful machine. 

The number of operations per step should be typical of explicit finite 
difference operators such as MacCormack^^^ 

The entropy production shown in Figure 1 illustrates that this numerical 
approximation to the variational statement does not guarantee monotonically 
increasing entropy and therefore does not absolutely guarantee stability. It 
is possible, however, that with a small enough time step or a higher order in 
time solution that monotonicity could be achieved (for a given problem). 
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Figure 4. - Jet Interaction 


4.0 CONCLUSIONS 


The foregoing discussion lends credence to the hypothesis that a 
variational principle for fluid mechanics does exist and that it can be 
utilized to create a powerful new approach to numerical solutions. The 
procedure developed here is applicable to problems of a complex yet practical 
nature which would be extremely difficult to solve by any other procedure. 

To summarize the developments discussed earlier, an analog has been 
created which is: 

o absolutely conservative 
o absolutely stable 

o explicitly adds no artificial viscosity 

o dynamically differences based on local conditions 

o can provide inexpensive and quick response for 
engineering and/or parametric solutions to whatever 
level of sophistication is deemed reasonable by the 
user. 

There is, of course, much more to be done to further verify the 
technique, investigate higher order solutions, and to extend it to general 
quadrilaterals, axisymmetric flows, three-dimensional flows, viscous flows, 
and reacting flows. 
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